library(lwdataexplorer)

zoo <- getZooscanData("2021-01-01", "2022-01-01")

head(zoo)

library(sf)

class(zoo)
#> [1] "data.frame"

zoo_sf <- st_as_sf(zoo, coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE)

class(zoo_sf)
#> [1] "sf"         "data.frame"

library(mapview)
mapview(zoo_sf)

# Access ETN with lwdataexplorer
etn <- getEtnData()

library(worrms)

taxa <- wm_records_taxamatch("Caretta caretta")

taxa[[1]]

library(mregions2)

# Find all region names starting with "Belgian"
belgian_gaz <- mr_gaz_records_by_names("Belgian")
belgian_gaz

# Get only Belgian Territorial Sea
library(mapview)
territorial_sea <- mr_gaz_record(49010)
mapview(territorial_sea)

library(sf)

nrow(zoo_sf)
#> [1] 437

zoo_sf_belgium <- st_intersection(zoo_sf, territorial_sea)

nrow(zoo_sf_belgium)
#> [1] 251

mapview(list(zoo_sf, territorial_sea))

mapview(list(zoo_sf_belgium, territorial_sea))

# API libraries
library(lwdataexplorer)
library(worrms)
library(mregions2)

## Helper libraries
library(mapview)
library(sf)


#---------------------------------------------------------------------------------------------
# Exercise 1.1.: 
# - Get fish tracking data from the European Tracking Network of the year 2015
#
# Hint: use `?getEtnData()` to find more information
#---------------------------------------------------------------------------------------------
etn <- getEtnData(
  startdate = "2015-01-01",
  stopdate = "2016-01-01",
  action = "Time bins",
  by = "1 day",
  networks = "All",
  projects = "All"
)

# Inspect
View(etn)


#---------------------------------------------------------------------------------------------
# Exercise 1.2.: 
# - Extract a list of species of the ETN data and 
# - Find more taxonomic information using the WoRMS Taxon Match from the worrms R package.
# - Left join the taxonomic info to the ETN dataset on the scientific name
#
# Hint: use `dplyr::distinct()` or `base::unique()` to get unique values
# Hint: worrms functions start all like `wm_*`
# Hint: use `do.call(rbind, my_data_frame)` to turn a list into a data frame
#---------------------------------------------------------------------------------------------

# Unique list of species
species_list <- unique(etn$scientific_name)

# Taxon match
species_matched <- wm_records_taxamatch(species_list)
species_matched <- do.call(rbind, species_matched)

# Left join
etn <- merge(etn, species_matched, by.x = "scientific_name", by.y = "scientificname")

# Inspect
View(etn)


#---------------------------------------------------------------------------------------------
# Exercise 1.3.: 
# - Use the mregions2 package to find and get the record for the Belgian Part of the North Sea.
# - Do a geospatial intersection between the latitude and longitude fields in the ETN 
#   data and the North Sea from Marine Regions. 
#
# Hint: you must transform the ETN data from a data.frame to a simple feature `sf` class. 
#       Use the function `st_as_sf()` with the `coords` and `crs` arguments
# Hint: all functions from the sf package start as `sf_*()`
# For example: `st_intersection()`
#---------------------------------------------------------------------------------------------

# Find the North Sea
mr_belgian <- mr_gaz_records_by_names("Belgian")
View(mr_belgian)

bpns <- mr_gaz_record(3293)

# Transform data frame into simple feature object
etn <- st_as_sf(etn, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

# Inspect both
mapview(list(etn, bpns))

# Perform the intersection
etn <- st_intersection(etn, bpns)

# Inspect again
mapview(list(etn, bpns))


#---------------------------------------------------------------------------------------------
# Bonus Exercise 1.4.: 
# - Save the ETN data into disk
#
#---------------------------------------------------------------------------------------------

# To turn into a data.frame again
st_drop_geometry(etn)

# Save as shapefile or csv
st_write(etn, "./data/etn.shp")
write.csv(st_drop_geometry(etn), "./data/etn.csv")

library(eurobis)
df <- eurobis_occurrences_basic(aphiaid = 126436, mrgid = 3293)
head(df)

mapview(df)

# base R
df$datecollected <- as.Date(df$datecollected)

# dplyr
library(dplyr)
df <- df %>%
  mutate(datecollected = as.Date(datecollected))

# API library
library(eurobis)

# Helper library
library(dplyr)

#---------------------------------------------------------------------------------------------
# Exercise 2.1.: 
#  - Use the eurobis R package to get EurOBIS data for the same species, time and location. 
#
# Hint: more info with ?eurobis_occurrences_basic()
#---------------------------------------------------------------------------------------------
mrgid = 3293
list_aphia = unique(etn$valid_AphiaID)

basic <- eurobis_occurrences_basic(
  mrgid = mrgid,
  aphiaid = species_matched$AphiaID,
  start_date = "2015-01-01",
  end_date = "2016-01-01"
)

#---------------------------------------------------------------------------------------------
# Exercise 2.2.: 
#  - Include more regions in your query! 
#    - Select your desired area with`eurobis_map_draw()` and pass to the argument `geometry` 
#    - Find more standardized marine regions with the family of functions `eurobis_map_regions_*()`
#  - Find data for the species of your interest by passing a scientific name.
#  - Find data for the functional group of your interest.
#---------------------------------------------------------------------------------------------
my_polygon <- eurobis_map_draw()

eurobis_map_regions_ecoregions()
eurobis_map_regions_eez()
eurobis_map_regions_iho()
eurobis_map_regions_eez_iho()
eurobis_map_regions_reportingareas()

basic <- eurobis_occurrences_basic(
  geometry = my_polygon,
  mrgid = c(2350),
  functional_groups = c("birds", "mammals")
)


#---------------------------------------------------------------------------------------------
# Exercise 2.3.:
#   - Wrangle the ETN dataset to fit the column names of eurobis.
#   - Select the following columns: datecollected, decimallongitude, decimallatitude, 
#     scientificname, aphiaid, scientificnameaccepted
#   - As we are not considering the individual fishes tracked but only presences, there will be 
#     duplicates. Get unique combination of values
#   - Bind the rows of both datasets
#
# Hint: the aphiaID from EurOBIS comes with a full URL. Use `gsub()` or `stringr::str_replace()` 
#       together with `as.numeric()` to get only the AphiaID 
# Hint: we recommend to use `dplyr::transmute()`, `unique.data.frame()` and `dplyr::bind_rows()`
# Hint: EurOBIS follows the Darwin Core standard for naming the columns: https://dwc.tdwg.org/terms/
#---------------------------------------------------------------------------------------------
colnames(etn)
colnames(basic)

basic <- basic %>% select(-gml_id, -id, -datasetid, -coordinateuncertaintyinmeters
             ) %>% mutate(datecollected = as.character(datecollected),
                          aphiaid = gsub("http://marinespecies.org/aphia.php?p=taxdetails&id=", "", aphiaid, fixed = TRUE)
             ) %>% mutate(aphiaid = as.integer(aphiaid))

df <- etn %>% 
  transmute(
    datecollected = time,
    decimallongitude = longitude,
    decimallatitude = latitude,
    scientificname = scientific_name, 
    aphiaid = AphiaID,
    scientificnameaccepted = valid_name
  ) %>% arrange(
    scientificname, datecollected
  ) %>%
  unique.data.frame(
  ) %>% 
  bind_rows(basic)


#---------------------------------------------------------------------------------------------
# Bonus Exercise 2.4.: 
#  - Repeat the query adding all the information available in EurOBIS
#  - Get all the distinct values of the parameters_* coumns
#  - Filter only records with count of individuals
#
# Hint: use `eurobis_occurrences_full_and_parameters()`
# Hint: more information here: https://www.emodnet-biology.eu/emodnet-data-format
#---------------------------------------------------------------------------------------------
full_emof <- eurobis_occurrences_full_and_parameters(
  mrgid = mrgid,
  aphiaid = species_matched$AphiaID,
  start_date = "2015-01-01",
  end_date = "2016-01-01"
)

# View all parameters available
full_emof %>% 
  st_drop_geometry() %>%
  select(parameter, 
         parameter_bodcterm, 
         parameter_bodcterm_definition,
         parameter_measurementtypeid) %>%
  distinct() %>% 
  View()


# Query only Counts
full_emof <- full_emof %>%
  filter(parameter_measurementtypeid == "http://vocab.nerc.ac.uk/collection/P01/current/OCOUNT01/")

library(EMODnetWFS)

emodnet_wfs()

# First create a WFS client
client <- emodnet_init_wfs_client("seabed_habitats_general_datasets_and_products")

class(client)
#> [1] "WFSClient"         "OWSClient"         "OGCAbstractObject" "R6"

# And ask for all the layers available
info <- emodnet_get_wfs_info(client)
info

# Download sandbank layer
emodnet_layer <- emodnet_get_layers(client, layers = "art17_hab_1110", 
                            crs = 4326, reduce_layers = TRUE)
class(emodnet_layer)
#> [1] "sf"         "data.frame"

head(emodnet_layer)

# Some editing needed to be able to plot
emodnet_layer <- emodnet_layer %>% 
  st_cast(to = "GEOMETRYCOLLECTION") %>% 
  st_collection_extract(type = "POLYGON")

mapview(emodnet_layer)

library(sdmpredictors)
library(raster)

# Find all layers available from the latest version
list_layers("Bio-ORACLE", version = 2.2)
list_layers_future("Bio-ORACLE", version = 2.2)

# Pass layer code for max sea surface temperature
load_layers("BO22_RCP26_2050_chlomean_ss")

list_layers("Bio-ORACLE", version = 2.2)